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SUMMARY 

The eigenvalue spectrum associated with a linear finite-difference approximation 
plays a crucial role in the stability analysis and in the actual computational perfor- 
mance of the discrete approximation. We investigate the eigenvalue spectrum associ- 
ated with the Lax-Wendroff scheme applied to a model hyperbolic equation. For an 
initial-boundary-value problem (IB VP) on a finite domain, the eigenvalue or normal 
mode analysis is analytically intractable. A study of auxiliary problems (Dirichlet and 
quarter-plane) leads to asymptotic estimates of the eigenvalue spectrum and to an 
identification of individual modes as either benign or unstable. The asymptotic analysis 
establishes an intuitive as well as quantitative connection between the algebraic tests 
in the theory of Gustafsson, Kreiss, and Sundstrom and Lax-Richtmyer L 2 stability on 
a finite domain. 

1. INTRODUCTION 


A classical method for carrying out a stability analysis of a discrete hyperbolic IB VP 
is the normal mode analysis of Gustafsson, Kreiss, and Sundstrom (GKS) [1], The GKS 
theory avoids the analytical intractability of the finite-domain normal mode analysis 
by analyzing related quarter-plane problems. On the other hand, when one performs 
numerical experiments to verify stability and/or accuracy predictions, the computations 
are on a finite domain and one typically uses the discrete L 2 norm and not the GKS 
norm used to prove stability. Thus in practice, we have the dichotomy of analyzing 
quarter-plane problems with GKS norms and computing on finite domains with L 2 
norms. 

The goal of this paper is twofold. First we present asymptotic limits for the normal 
modes of the discrete (Lax-Wendroff) IB VP on a finite domain. These limits lead to a 
delineation of the normal modes of the finite-domain problem into three classes. Next 
we use the asymptotic estimates to make a direct algebraic connection between the 
normal modes of the GKS quarter- plane analysis and the classes of normal modes of 
the finite-domain problem. This leads to an interpretation of (unstable) GKS modes 
which is readily understandable in terms of the Lax-Richtmyer stability in the L 2 norm. 
In this paper we give only a brief outline of our analysis. A detailed exposition is given 
in [2]. 

2. IBVP FOR A MODEL HYPERBOLIC EQUATION 


We consider the scalar hyperbolic equation 


du du 
dt dx ’ 


0 < x < L, 


t> 0 


( 2 . 1 ) 
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wljere u = u(x,t ) and c is a real constant. For a well-posed IBVP on a finite domain 
one must specify initial data, ti(x,0) = /(x), 0 < x < Z-, and an analytical boundary 
condition at x = L, 

u(L,t) = g(t) for c > 0. (2.2) 

3 7 . A PROTOTYPE FINITE-DIFFERENCE APPROXIMATION 

To obtain a difference approximation of the model equation (2.1) a mesh is introduced 
in (x, t) space with increments Ax and At and indexing defined by x = jiAx and 
t = uAt. The spati.al domain 0 < x < T is divided into J equally spaced increments, 
i.e., J Ax — L. As a prototype (explicit) finite-difference approximation for the model 
equation (2.1), we consider the Lax-Wencfroff scheme 

u r' = + § (“?+. - + yK+. - K + <*?_!> (3.1) 

where v, = cAtj Ax is defined to be the Courant number. In our analysis the an- 
alytical boundary conditipii (2.2) for the difference approximation is assumed to be 
homogeneous, i.e., 

u 1 } = 0. (3.2) 

If we apply the Lax-Wendroff scheme at the outflow boundary (j = 0), the com- 
putational stencil protrudes one point to the left of the boundary. It is clear that an 
additional numerical boundary scheme (NBS) is required to calculate i.e., the 

solution pn the outflow boundary at time level n + 1. As a prototype NBS we choose 
the spatially one-sided scheme: 

< +1 = + v[—au% + (1 + 2 a)u™ — (1 + a)uQ ] (3*3) 

where a is a (real) parameter. If a = 0, then (3.3) is simply 

< +1 = < + v(ul - <). (3.4) 

The bax-Wendroff scheme (3.1) together wi.tb.the analytical boundary condition (3.2) 
and the NBS (3.3) is.called a discrete IBVP. For our purposes it. is convenient to rewrite 
tfie NBS (3.3) as an equivalent space extrapolation formula [2]: 

h(E)uly= 0 where h{E) = (E - l) 2 [2aE - (1 - u)). (3.5) 

The. shift operator. E is defined by Euj = uj + 1 and h(E) is a polynomial in E- 

4.. LAX-RICHTMYER stability of a discrete ivp OR IBVP 

For the. stability analysis of a discrete initial value problem (IVP) or IBVP with 
requisite homogeneous boundary conditions, it is appropriate to write the discrete 
apprpximartion in vector-matrix form: 

u n+1 = Cu n . (4.1) 
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The Lax-Wendroff scheme (3.1) with the analytical boundary condition (3.2) and 
NBS (3.3) can be written in vector- matrix form (4.1) where 




m s V w 

u ? 


p q r 0 

• 

, C = 

« * * 

• 


. 

U J - 2 


0 p q r 

- U J - 1 - 


1 


(4.2a, b) 


and 

p = —i/(l — v)/ 2, q = 1 — v 2 , r = i/(l + i/)/2 (4.3a) 

s = 1 - i/(l + a), v = i/(l + 2a), w = -va. (4.3b) 

Here the matrix size is J x J. 

The determination of the eigensolutions of the first-order system (4.1) is sometimes 
called the normal mode analysis. If C has a complete set of eigenvectors, then the 
general solution of (4.1) can be written as 


J - 1 

u n = 

l- 0 


(4.4) 


where zi and <f>( denote the £th eigenvalue and eigenvector of the matrix C and the a^’s 
are complex constants determined from a specified initial vector u°. Thus the eigenso- 
lutions are the normal modes and it is obvious from (4.4) that they act independently. 

A discrete IVP or IB VP represented by (4.1) is Lax-Richtmyer stable if there exists 
a constant K > 1 such that for any initial condition u° 


||u"|| < A||u°|| (4.5) 

for all n > 0, 0 < nAt < T with T fixed and At/ Ax fixed. A necessary condition for 
Lax-Richtmyer stability is that there exists a nonnegative constant w such that 

P(C) < 1 + j (4.6) 

for all n > 0, 0 < nAt < T with T fixed and At/ Ax fixed. Here />(C) denotes the 
spectral radius of C. Inequality (4.6) is referred to as the spectral radius condition. 

5. EIGENVALUE SPECTRUM OF A DISCRETE IVP 

A necessary condition for the stability of a discrete IBVP is the stability of the 
corresponding pure IVP or Cauchy problem. In this section we review the eigenvalue 
spectrum of the IVP for the Lax-Wendroff scheme. The solution of the discrete IVP is 
assumed to be spatially periodic with period L — J Ax, and hence 


U 1 = u" 




(5.1) 
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Consequently, the Lax-Wendroff scheme can be written in vector-matrix form (4.1) 
where C is a J x J circulant matrix and the eigenvalues Z( are given analytically by 

z t — 1 - 2i/ 2 sin 2 (0*/2) + ivsinOt, l = 0, 1, • • • , J - 1 (5.2) 

where = 2£ir j J and v = cAt/ A® is the Courant number. The eigenvalue locus given 
by (5.2) is an ellipse in the complex z-plane. (The eigenvalue locus is defined to be a 
curve through the eigenvalues in the complex z-plane.) The eigenvalue loci for v = 0.5 
and v — 1.1 are shown in Fig. 5.1. 



Fig. 5.1. Eigenvalue locus for periodic and Dirichlet boundary conditions. 


The ellipse is contained within the unit circle for \u\ < 1 (Fig. 5.1a) and if \u\ = 1, 
the locus is the unit circle. If \u\ > 1, the ellipse contains the unit circle (Fig. 5.1b). 
Since a circulant matrix is normal, the spectral radius condition (4.6) is necessary and 
sufficient for stability. Consequently, as is well-known, the Lax-Wendroff scheme is 
stable for the IVP (i.e., Cauchy stable) in the L 2 norm for It'] < 1. 

8. EIGENVALUE SPECTRUM OF A DISCRETE IBVP 


Before we present asymptotic limits for the eigenvalues of the IBVP matrix (4.2b) 
it is advantageous to examine the eigenvalue spectrum by computing the eigenvalues 
numerically. As an example if we choose a = 0 in (3.3) we obtain the NBS (3.4). A 
sketch of the eigenvalue locus is shown by the dashed curve in Fig. 6.1a for v = 0.5. 




(a) z-plane (a = 0) (b) z-plane (a = -1) 

Fig. 6.1. Dashed curve is eigenvalue locus of the matrix (4.2b) for i/=0.5. 
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If we increase the number of spatial increments J one finds that the eigenvalue locus 
approaches the solid vertical line for increasing J . 

As a second numerical example, we consider a = —1 in (3.3). This NBS leads to an 
unstable discrete IBVP. The eigenvalue locus is shown in Fig. 6.1b. The eigenvalue lo- 
cus is qualitatively similar to the locus shown in Fig. 6.1a except for a single eigenvalue 
outside the unit circle shown by the solid symbol in the figure. To plotting accuracy 
the eigenvalue indicated by the solid symbol appears fixed for increasing values of J . 
The eigenvalue locus indicated by the dashed curve approaches the solid vertical line 
for increasing values of J. If one does a GKS quarter-plane analysis, one finds that the 
lone eigenvalue outside the unit circle is a GKS eigenvalue. 

7. AN AUXILIARY DIRICHLET PROBLEM 

The GKS stability analysis involves three auxiliary problems: the Cauchy problem 
and the left- and right-quarter plane problems. In this section we consider a fourth 
auxiliary problem which we call the Dirichlet problem on a finite domain. The impor- 
tance of the Dirichlet problem accrues from the fact that, with the exception of isolated 
eigenvalues detected by the GKS quarter-plane analysis, the eigenvalue spectrum of a 
discrete IBVP is a perturbation of the eigenvalue spectrum of the auxiliary Dirichlet 
problem. 

The auxiliary Dirichlet problem is constructed by equating to zero any grid function 
value Uj which is required by the interior difference approximation but falls outside 
the computational domain 0 < x < L. Hence the auxiliary Dirichlet problem for the 
Lax-Wendroff scheme can be written in the vector-matrix form (4.1) where the matrix 
operator C is a J x J tridiagonal matrix with elements p, g,r. The eigenvalues of C 
can be determined analytically: 

zt = (1 — v 2 ) + iv\J\ — v 2 cos 0/, £ = 1, 2, • ■ • , J (7.1) 

where = h r/(J + 1). If \u\ < 1, the eigenvalues are complex and the eigenvalue locus 

is a vertical line centered at the point (1 — z^ 2 , 0) in the complex z-plane. If \u\ = 1 all 
the eigenvalues degenerate to the single point zt — 0. For \v\ > 1 the eigenvalues are 
real. The eigenvalue spectra of the pure IVP (periodic boundary conditions) and the 
auxiliary Dirichlet problem are compared in Fig. 5.1. By some elementary calculations 
[2] one finds the rather remarkable result that the eigenvalue locus of the auxiliary 
Dirichlet problem is simply a straight line segment joining the foci of the ellipse which 
is the eigenvalue locus of the IVP. 

8. NORMAL MODE ANALYSIS 

In this section the relevant formulas for the normal mode analysis of the finite-domain 
problem and the quarter-plane problem are summarized. A detailed analysis is given 
in [2]. 

8.1 Finite-Domain Normal-Mode Analysis — Summary. An eigensolution or 
normal mode of the finite-domain IBVP is determined by looking for a solution of the 
form 

u] = ( 8 . 1 ) 
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which satisfies the Lax-Wendroff scheme (3.1) with the analytical boundary condition 
(3.2) and the NBS (3.3) written as an extrapolation formula (3.5). 

The eigenvalue z is given by 


„ v , lx v 2 , 1. 

Z=z1 + 2 {k ~k ) + Y (k ~ 2+ k ) (8<2) 

and the components <f>j of the eigenvector <f> are 

4>j = aW ~ (-« 2 /C) J (-C/^] (8.3) 

where « = y/(k and k is a root of the characteristic equation 

M\/C«)-(-« 2 ) J+1 M-v / CM) = o. (8.4) 

The parameter (, which is positive for — 1 < v < 1, is defined by 


c = 


1 - V 

T+V 


(8.5) 


The coefficient a on the right-hand side of (8.3) is an arbitrary constant. The polynomial 
h(^k) depends solely on the NBS (3.5), i.e., h(y/£k) is the polynomial associated 
with the NBS written as an extrapolation formula. If one could solve for the roots of 
the characteristic equation (8.4), then the eigenvalues z and the eigenvectors <f> would 
follow directly from (8.2) and (8.3). The normal mode analysis on a finite domain 
is, in general, analytically intractable because one cannot solve for the roots of the 
characteristic equation (8.4). 

For the auxiliary Dirichlet problem, the polynomial h(y/£k) is unity and (8.4) reduces 
to 

(-« 2 ) J+1 = 1 . ( 8 . 6 ) 

One can solve (8.6) by using the roots of unity formula and the normal modes can be 
found analytically. 

8.2 Quarter-Plane Normal-Mode Analysis - Summary. For the right quarter- 
plane problem one also looks for a solution of the form (8.1) which satisfies the Lax- 
Wendroff scheme (3.1) and the NBS (3.5). The details of the GKS normal mode analysis 
are given in [2, Appendix A], The eigenvalue 2 is given by (8.2) and the components 
4>j of the eigenvector <j> are 

4>j = (8.7) 

where k = y/^k and k is a root of the quarter-plane characteristic equation 

h(V(k) = 0. (8.8) 
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Fig. 8.1. GKS-stability region for Lax-Wendroff (3.1) with NBS (3.3). 

The (a,^) parameter space for which the discrete IBVP is GKS stable is shown by the 
cross-hatched region of Fig. 8.1. 

9. ASYMPTOTIC ROOTS OF CHARACTERISTIC EQUATION 

There exists a small class of discrete IBVP’s which are sometimes called border- 
line cases. Borderline cases are unstable according to the GKS theory but they may 
be Lax-Richtmyer stable or unstable in the L 2 norm on a finite domain. Borderline 
approximations can be characterized by the presence of a stationary mode for the finite- 
domain problem. A stationary mode is defined to have the property that a k root of the 
characteristic equation (8.4) is independent of J, i.e., k remains fixed in the complex 
plane as J increases. It is important to note that the detection of a stationary mode 
requires no asymptotic (large J ) analysis because there is no J dependence. 

Although there can be stationary modes which are independent of J, almost all roots 
of the characteristic equation (8.4) depend upon J and we write k = k(J). One can 
show that there is no loss in generality in assuming |/c| < 1 and we write 

|/c| = |k(J)| = 1 — e, 0 < e(J) < 1 (9.1) 


where 

either e(J) >6> 0 or e(J) — ► 0 as /-too. (9.2a, b) 

In particular we are interested in the conditions under which the characteristic equa- 
tion (8.4) reduces to the quarter-plane characteristic equation (8.8) in the limit J — > 00 . 
Obviously this depends on the asymptotic behavior of |/c( J )| J . There are only two pos- 
sibilities, either 

/: lim \k{J)\ J = constant > 0, or II : lim |/c(/)| J = 0. (9.3a,b) 

In case II it is obvious that (8.4) reduces to (8.8) as J — + 00 . In case I one can show 
that the roots of (8.4) asymptotically approach the roots of the auxiliary Dirichlet 
problem (8.6) as J — + 00 . The details of the asymptotic estimates are in [2]. 
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10. CLASSIFICATION OF NORMAL MODES 

Tft€ normal modes of a discrete IBVP can be divided into three classes according to 
their asymptotic behavior as / — » oo. We associate the following nomenclature with 
the three classes: 

-f. Dirichlet- h&e inodes 

II. Quarter-plane- life modes (10-1) 

HI- Stationary modes. 

The adjective like is used to imply that one can identify for finite / a particular mode 
tjmt becomes either a Dirichlet mode or quarter-plane mode in the limit J — > oo. 

Fqr a given difference approximation almost all the normal modes are in class I. 
These modes have a generic eigenvalue spectrum which is easy to describe; to wit, for 
a fPye.u Cpo.rant number y and finite / the spectrum is. simply a perturbation of the 
Diriclilet spectrum. A typical generic ease is shown in Fig. 6.1a. The solid vertical 
line is the Djrichlet locus. The dotted line slightly to. the left is the eigenvalue locus of 
the cla^s I inodes of the difference approximation to the IBVP. One can show [2] that 
the dotted eigenvalue locus approaches the Dirichlet locus at least as fast as 0(1//). 
Hence as the mesh is refined, i.e., / > oo, the dotted, locus collapses onto the Dirichlet 

locus- Hence the terminology Dirichlet- like modes. 

The class I modes are always benign in the sense that they do not introduce unstable 
mode.s into a difference approximation which is Cauchy stable. Only modes q£ class II 
aod III can, introduce unstable modes into a difference approximation, which, is Cauchy 
If modes of class II apd III exist, they are created by the NBJS. 

Tl^e modes ip class 1/ are related to the G.KS stability theory in the sense that they 
become quarter-plane modes as / — > oo, i.e., as the the mesh is refined. Finally, the 
stationary modes which constitute cfass^ III arq common to. both the finite-domain 
problem and the quarter-plane problems. 

11. SKETCHES OF ROOT AND EIGENVALUE DISTRIBUTIONS 

In t)ds section we give a pictorial description of. the roots of the characteristic equation 
(8.4) and, the corresponding eigenvalue spectra, associated with each of the three classes 
(t0»l). The k roots plotted in the examples were computed numerically from the 

c ;I iara £teristic equation (8.4) and the corresponding eigenvalues were computed using 

' ' 

The r . roots for the characteristic equation (8.6) of the auxiliary Dirichlet problem. are 
Pitted ip Fig. 11.1^ % / = 19, The corresponding eigenvalue, locus is the vertical 
lipe .sl^own in Fig. 11.1b. In this figure and tile figures ,to follow we plpt the eigenvalue 
locus^rather, than Individual eigenvalues because of the small size of the figures. 

Ill Qirjchlet- like modes. The roots, of the characteristic, equation (8.4) which 
correspond tp modes in class I have a generic, root locus in the complex /c-plane which 
if simply a perturbation (inside, the unit circle) of the Dirichlet. root locus which is on 
t^ui^l, circle. As an .example we consider a stable case from.the shaded region of Fig. 
8.1 by choosing parameter values a = 0 apd .^ — 0.5. The roots of (8.4) for / = 19 are 
plotted in Fig. 11.2a. 
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(a) «-plane 

Fig. 11.1. Roots of characteristic polynomial 



(8.6) and eigenvalue locus for v = 0.5. 



(a) /c-plane (a = 0) 

Fig. 11.2. Roots of characteristic polynomial 



(8.4) and eigenvalue locus for v = 0.5. 


The eigenvalue locus is indicated by the dashed curve of Fig. 11.2b. The solid vertical 
line is the eigenvalue locus of the auxiliary Dirichlet problem. As J increases the root 
locus approaches the unit circle and the eigenvalue locus moves toward the Dirichlet 
locus at least as fast as 0(1/ J). 

11.2 Quarter- plane- hie modes. The previous example had only class I modes, 
i.e., Dirichlet - like modes, while the examples of this section have both class I and class 
II modes. Modes in class II are related to the GKS theory in the sense that they 
become quarter- plane modes as J — > oo. In addition, there are only a few modes in 
this class and the maximum number is known exactly. 

The following two examples are unstable discrete IBVP’s. Even though a discrete 
IBVP is unstable, there is no dramatic change in the eigenvalue locus in the sense that 
it remains a perturbation of the Dirichlet eigenvalue locus but with the addition of one 
or two eigenvalues near or strictly outside the unit circle. These additional eigenvalues 
correspond to GKS eigenvalues or generalized eigenvalues in the limit J — > oo. 

11.2a Unstable quarter-plane- like mode — GKS eigenvalue. For the first 
example we choose a = — 1 and v = 0.5 from the unshaded region of Fig. 8.1. The 
roots of (8.4) are shown in Fig. 11.3a for J — 19. 

From the figure it is apparent that there is a single isolated root indicated by the solid 
symbol in the figure. The corresponding eigenvalue is indicated by the solid symbol in 
Fig. 11.3b. This single eigenvalue remains strictly outside the unit circle as J — + oo 
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(a) /e-plane (a = —1) (b) z-plane (a = — 1) 


Fig. 11.3. Roots of characteristic polynomial (8.4) and eigenvalue locus for v = 0.5. 

and consequently the approximation is unstable (GKS eigenvalue). 

11.2b Unstable quarter- plane- hie mode - GKS generalized eigenvalue. One 

should expect the discrete IBVP to be unstable for u < 0 for otherwise one would have 
a stable approximation for an ill-posed IBVP. As an example we choose a — 0 and 
v = —0.5 from the unshaded region of of Fig. 8.1. The roots of (8.4) are depicted in 
Fig. 11,4a for J = 19. 



(a) /c-plane (a = 0) (b) z-plane (a = 0) 

Fig. 11.4. Roots of characteristic polynomial (8.4) and eigenvalue locus for v — —0.5. 

There are two (nearly coincident) isolated roots as shown by the solid symbol of the 
figure. The corresponding eigenvalues are z ^ 1 in Fig. 11.4b. For finite J, one of these 
eigenvalues is slightly inside the unit circle and the other is slightly outside. In this 
example the origin of the instability is rather subtle. The instability is not due to a 
violation of the spectral radius condition (4.6) but rather is due to the introduction of 
a solution (proportional to the eigenvector) whose norm cannot be uniformly bounded 
by the norm of the initial data as the mesh is refined, i.e., there is algebraic growth. In 
the nomenclature of the GKS theory, there is a generalized eigenvalue. 

11.3 Stationary mode. As a final example we consider a borderline case. For pa- 
rameter values we pick a = —0.75 and u = 0.5 which is a point on the left boundary 
between stability and instability in Fig. 8.1. The roots of (8.4) are shown in Fig. 11.5a 
for J = 19. 
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Fig. 11.5. Roots of characteristic polynomial (8.4) and eigenvalue locus for v = 0.5. 


There is one isolated root shown by the solid symbol of the figure. The corresponding 
eigenvalue is z = 1. This eigenvalue of unity is independent of J. Since z n = 1, the 
stability or instability of the difference approximation devolves to the behavior of the 
corresponding eigenvector as the mesh is refined. This approximation happens to be 
Lax-Richtmyer stable but GKS unstable. The details of this example are worked out 
in [3]. 

12. CONCLUSIONS 

We have investigated the eigenvalue spectrum for the Lax-Wendroff scheme applied 
to a model hyperbolic IBVP. For the discrete IBVP on a finite domain, the normal 
mode analysis is analytically intractable even for the simple prototype difference ap- 
proximation. On the basis of an asymptotic normal mode analysis (large J ), we have 
classified the normal modes of the finite-domain problem into three classes. The result- 
ing classification leads to a simple description of the asymptotic eigenvalue distribution. 
For a given Courant number, the spectrum is simply a perturbation of the spectrum 
of the auxiliary Dirichlet problem plus whatever eigenvalues are present in the GKS 
normal mode analysis for the related quarter-plane problems. 

Almost all the modes are Dirichlet - like modes and are benign in the sense that they do 
not introduce unstable modes into a difference approximation which is Cauchy stable. 
Only quarter-plane- like modes and stationary modes can introduce unstable modes into 
an approximation which is Cauchy stable. Consequently, if an instability exists it is 
caused by the NBS and is detected in the GKS analysis. 
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